Strong polymer-turbulence interactions in viscoelastic turbulent channel flow 
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This paper is focused on the fundamental mechanism(s) of viscoelastic turbulence that lead to 
polymer induced turbulent drag reduction phenomenon. A great challenge in this problem is the 
computation of viscoelastic turbulent flows, since the understanding of polymer physics is restricted 
to mechanical models. An effective state-of-the-art numerical method to solve the governing equa- 
tion for polymers modelled as non-linear springs, without using any artificial assumptions as usual, 
was implemented here for the first time on a three-dimensional channel flow geometry. The capa- 
bility of this algorithm to capture the strong polymer-turbulence dynamical interactions is depicted 
on the results, which are much closer qualitatively to experimental observations. This allowed a 
more detailed study of the polymer-turbulence interactions, which yields an enhanced picture on a 
mechanism resulting from the polymer-turbulence energy transfers. 



I. INTRODUCTION 

A few parts per million by weight polymer in a wall- 
bounded turbulent flow are enough to reduce the force 
necessary to drive the flow through a channel by a fac- 
tor of up to 70%, as was discovered by Toms [l| while 
performing experiments on the degradation of polymers. 
Turbulence is a multiscalc phenomenon with a vast spec- 
trum of spatial scales and therefore a very large num- 
ber of degrees of freedom. Due to the fact that even 
the maximum polymer molecule end-to-end distance Lp 
is much less than the Kolmogorov viscous scale 77, one 
might anticipate that the small size polymers can only 
affect sub-Kolmogorov scale processes and that scales of 
length £ > i] would remain unaffected. Surprisingly, the 
dynamics of the small polymer chains arc able to funda- 
mentally modify the large scale structures and statistics, 
as observed in the drag reduction (DR) phenomenon Q . 

Polymer drag reduction in wall-bounded turbulent 
flows induces higher mean velocities, implying deviations 
from the classical phenomenology of hydrodynamic wall- 
bounded turbulence and hence from the von Karman law 
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where '+' denotes normalisation with the friction veloc- 
ity u r and the viscous length scale S v = vju T with v 
being the fluid's kinematic viscosity. Moreover, k in Eq. 
(QJ is the von Karman coefficient [U, Q, which is usu- 
ally considered to be a constant taking the value 0.41 
and B ~ 5.2 is the intersept constant. The detailed ex- 
perimental work by Warholic et al. [BJ distinguished be- 
tween polymer induced drag reduced flows at low drag re- 
duction (LDR) and high drag reduction (HDR) regimes, 
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based on the statistical trends of the turbulent velocity 
field. When |DR| < 40% (LDR), the mean velocity pro- 
file is a log-law parallel to the von Karman law Eq. (UJ 
with a higher value of B, i.e. larger mean velocity. How- 
ever, for 40% < |DR| < 60% (HDR), the slope of the 
log-region increases until it reaches the empirical maxi- 
mum drag reduction (MDR) asymptote 
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where k^ 1 ~ 11.7 and the intercept constant B v ~ 17. 
This universal asymptotic profile was discovered experi- 
mentally in pipe flow by Virk et al. @, 0] and confirmed 
experimentally in channel flow by Warholic et al. [f|. 
Overall, the mean velocity profile is bounded between 
the von Karman law Eq. JT]) and the MDR law Eq. ©, 
the latter being independent of the Newtonian solvent, 
the polymer characteristics and the flow geometry. 

Polymer induced drag reduction has been known for 
more than sixty years and has attracted attention both 
from the fundamental and applied perspective. However, 
no generally accepted theory has been provided to ex- 
plain adequately the phenomenon. Such a theory should 
provide an explanation of the drag reduction onset, as 
well as the MDR law and its universality, which plays 
a significant fundamental role in understanding the phe- 
nomenon. Several theoretical concepts have been pro- 
posed but all have been subjected to criticism. The pro- 
posed theories fall mainly into two categories, that of 
viscous [1, HJ and that of elastic effects flol - [l2j . 

Recent progress in DNS of viscoelastic turbulence has 
begun to elucidate some of the dynamical interactions 
between polymers and turbulence, which are responsible 
for drag reduction. The aim of this study is to investigate 
the polymer dynamics, their influence on flow quantities 
and the phenomenology of drag reduction in the various 
regimes through DNS of viscoelastic turbulent channel 
flow using the finite extensible non-linear elastic model 
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with the Peterlin linearisation (FENE-P) [15| . the most 
widely used coarse-grained model in such studies. 

The paper is organised as follows. The necessary de- 
tails on the DNS of viscoelastic turbulent channel flow 
are provided in section [TTJ We analyse various viscoelas- 
tic turbulent statistics in section IIIII for all the drag re- 
duction regimes achieved in this study with a state-of- 
the-art numerical approach, which we have adapted to 
wall-bounded flows [l3[, aimed at capturing discontinu- 
ities in the polymer field. This approach is described in 
some detail in appendix [TJ Specifically, the effects of 
polymer extensibility and Reynolds number are briefly 
considered, whereas the statistics of mean velocity, fluc- 
tuating velocities and vorticitics are examined in depth 
demonstrating that our computations are qualitatively 
closer to experimental observations than previous numer- 
ical studies. Section ITVl presents the conformation tensor 
statistics and the scaling of polymer stress tensor com- 
ponents at the high Weissenberg number limit, which 
assists in a new asymptotic result for the shear stress 
balance (see section [V}. Finally, the polymer-turbulence 
interactions are studied in section IVT1 through the energy 
balance. A refined and extended picture of a conceptual 
model for drag reduction based on viscoelastic dissipa- 
tion is proposed in section IVIII before summing up our 
most important results (see section |VIII|) . 



II. DNS OF VISCOELASTIC TURBULENT 
CHANNEL FLOW 



The enormous number of degrees of freedom of each 
coil means that polymers are an extraordinarily complex 
system, whose dynamics depend on the conformations 
of the polymer molecules, i.e. orientation and degree 
of stretching of a coil. The study of detailed motions 
of this complex system and their relations to the non- 
equilibrium properties would be prohibitive. Only af- 
ter elimination of the fast relaxation processes of local 
motions in favour of stochastic noise, is it possible to 
study the dynamics of longer relaxation time scales [3] , 
such as the end-to-end conformation, that are responsible 
for many physical properties of polymers in fluids, such 
as viscoelastic turbulence and polymer drag reduction. 
Thus, coarse-grained mechanical models, such as bead- 
rod-spring models, are very crucial in DNS of viscoelastic 
turbulence. 

The computationally demanding Navicr-Stokes equa- 
tions in three-dimensions makes a Lagrangian approach 
for the polymer equally prohibitive and also limits poly- 
mer models to simple representations. A successful model 
for DNS studies of turbulent drag reduction is the FENE- 
P model in the Eulerian frame of reference, representing a 
conformation field of polymer macromolecules that have 
been modelled as non-linear bead-spring dumbbells [lfj. 
The standard approach to numerically solve the FENE-P 
model and its slight variations [TH, [l7| add an artificial 
diffusion term in the conformation field equation to avoid 



the loss of strict positive definiteness (SPD) of the con- 
formation tensor and subsequently numerical breakdown 
caused by the hyperbolic nature of the FENE-P model 
(see Eq. ©). 

Jin and Collins [lH stress the fact that much finer grid 
resolutions are required to fully resolve the polymer field 
than the velocity and pressure fields. Indeed, the hyper- 
bolicity of the FENE-P model admits near discontinuities 
in the conformation and polymer stress fields [l9| . Qual- 
itatively similar problems occur with shock waves and 
their full resolution in gas dynamic compressible flows, 
which is not practical using finer grids. In this case, 
high resolution numerical schemes such as slope-limiter 
and Godunov-type methods [2(| have proved successful 
at capturing the shock waves by accurately reproducing 
the Rankine-Hugoniot conditions across the discontinuity 
to ensure the correct propagation speed. 

Motivated by these schemes, Vaithianathan et al. [2l[ 
adapted the second-order hyperbolic solver by Kurganov 
and Tadmor [22| , which guarantees that a positive scalar 
remains positive over all space, to satisfy the SPD prop- 
erty of the conformation tensor and therefore avoid loss 
of evolution. Vaithianathan et al. [2l] | further demon- 
strated that this scheme dissipates less elastic energy 
than methods based on artificial diffusion, resulting in 
stronger polymer-turbulence interactions. Moreover, ac- 
cording to the most recent review on the subject 0], there 
are a lot of divergent and misleading results because of 
this artificial term introduced in the governing equations. 
For these reasons a modification of this shock-capturing 
scheme was developed in this present study to comply 
with non-periodic boundary conditions (see appendix fl~|) . 



A. Governing equations 

The dimensionless incompressible Navier- Stokes equa- 
tions for a viscoelastic fluid take the form 



V u = 



d t u + (u ■ V)tt = - Vp + — Atj + V • er 

Re c 



(3) 



where f3 = fi s / is the ratio of the solvent viscosity fi s 
to the total zero-shear-rate viscosity of the solution /i , 
Re c = U c S/u is the Reynolds number based on U c = |£/ft 
with Ub the bulk velocity of the flow kept constant in time 
and the channel's half-width 5. The extra force in Eq. 
([3]) arises due to polymers and the polymer stress tensor 
for the FENE-P dumbbells is defined by the Kramers 
expression 



1-/3 
Re c We c 



(f(trC)C - I) 



(4) 



where We c = t p U c /S with t p the polymer relaxation time 

L 2 —3 I — I, 

scale, f(trC) = L -i v _ trC is the Peterlin function [23[ and 
C = (QQ) is the conformation tensor, which is defined 
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as the dyadic product of the end-to-end vector Q of a 
dumbbell that specifies its configuration. The Peterlin 
function prevents the dumbbell to reach its maximum 
extensibility, i.e. trC < L 2 , since as trC — > L 2 the force 
required for further extension approached infinity. Note 
that C and L 2 are made dimensionless by the equilib- 
rium length scale y/k^T/H, where fee is the Boltzmann 
constant, T is the solution temperature and H is the 
Hookean spring constant and they have been normalised 
such that the equibrium condition is C eq = I. Then, the 
conformation tensor is governed by the FENE-P model 



d t C+{u-V)C-C-Vu-Vu T -C 



We, 



■(f(trC)C-I) 

(5) 

where the left hand side is the material derivative for 
a tensor field preserving its Galilean invariance and the 
right hand side represents deviation from the isotropic 
equilibrium due to Warner's finite extensible non-linear 
elastic spring- force law [24| . 

The clastic potential energy per unit volume E p stored 
by FENE-P dumbbells can now be specified by taking the 
integral of Warner's spring-force law over the end-to-end 
vector and after some algebra we obtain 

E P = \^^i L l ~ 3 ) H.f(trC)) + E po (6) 

where E po is a constant reference energy at equilibrium. 
After that, taking the time derivative of the elastic po- 
tential energy 



d t E p = 
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using the trace of Eq. ((5j), viz. 
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and similarly for the V-E p , we can derive the follow- 
ing balance equation for the elastic potential energy of 
FENE-P dumbbells 



dtE p + u ■ Vi? p = er • Vu 
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where E p is produced by er • Vw, dissipated by 
2Wc f{t r C)trcr and transported by u ■ 'VEp. 

B. Numerical parameters and procedures 

Incompressible viscoelastic turbulence in a channel was 
simulated in a rectangular geometry by numerically solv- 
ing the non-dimensional Eqs. (JS])-© in Cartesian co- 
ordinates. After obtaining the new update of the con- 
formation tensor from the FENE-P model using the 
method described in appendix fTl Eqs. Q-© ar e numer- 
ically integrated with a fractional step method using a 



second-order Adams-Bashworth/Trapczoidal scheme (see 
appendix l"2"|) . The fractional step method projects the 
velocity field to a divergence free velocity field and the 
Poisson pressure equation is solved in Fourier space with 
a staggered grid for the pressure field [lH . The staggered 
grid for the pressure was used for numerical stability 
purposes as well as the skew-symmetric implementation 
of the non-linear term in Eqs. ([3]). Spatial derivatives 
are estimated using sixth-order compact finite-difference 
schemes [26|. The grid stretching technique used in the 
inhomogeneous wall-normal direction maps an equally 
spaced co-ordinate in the computational space to a non- 
equally spaced co-ordinate in the physical space, in or- 
der to be able to use Fourier transforms [25|, (27| ■ Fur- 
ther details of our numerical method are provided in 
[l3j. Moreover, a validation of the algorithm just for the 
Navicr-Stokes equations for turbulent channel flow can be 
found in }25j , where this methodology was compared with 
spectral and second-order finite-difference schemes show- 
ing the necessity of spectral-like accuracy of the compact 
high-order schemes against second-order finite-differences 
in turbulence computations. 

To simulate incompressible channel flow turbulence we 
applied periodic boundary conditions for u = (u, v, w) in 
the x and z homogeneous directions and no-slip boundary 
conditions u = at the walls. The mean flow is in the x 
direction, i.e. (u) = ((u(y)), 0, 0), where ( } in this paper 
denotes averages in x, z spatial directions and time. The 
bulk velocity Ub in the x direction was kept constant 
for all computations at all times by adjusting the mean 
pressure gradient —d(p)/ dx at each time step. The choice 
of Ub in the computations for the Newtonian fluid is made 

" 2l for 
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based on Dean's formula Rc To ~ 0.119Rcy s [2 



a required Rc T 



where u TQ is the friction velocity 



for Newtonian fluid flow, i.e. j3 — 1 (see N cases in Table 

Q. 

The procedure used for the computation of the vis- 
coelastic turbulent channel flows of Table fl] is the follow- 
ing. First, DNS of the Newtonian fluid, i.e. (j = 1, were 
performed for the various Reynolds numbers until they 
reached a steady state. Then, the initial conditions for 
the viscoelastic DNS were these turbulent Newtonian ve- 
locity fields as well as the stationary analytical solution of 
the FENE-P model, given a steady unidirectional shear 
flow u = (U(y), 0,0), for the CV, tensor components 
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with 4> = cosh^ 1 (f n 2 + 1), o = and a_ v 
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-6(y-l) 7 assuming that U(y) = 0.75(1 - (y - l) 8 ) Vy E 
[0, 2] is a close approximation to the averaged velocity 
profile of a Newtonian fully developed turbulent chan- 
nel flow at moderate Reynolds numbers [30T] . Initially, 
the governing equations were integrated uncoupled, i.e. 
(5 = 1, until the conformation tensor achieved a station- 
ary state. From then on the fully coupled system of equa- 
tions, i.e. P ^ 1, was marched far in time, while u and 
C statistics were monitored for several successive time 
integrals until a fully developed steady state is reached, 
which satisfies the total shear stress balance across the 
channel, viz. 



/? d(u) 
Re c dy 




where (a 12 ) = r^w^ ( L 2l Cfcfc C12) is the mean poly- 



mer shear stress. Finally, after reaching a statistically 
steady state, statistics were collected for several decades 
of through-flow time scales L x /Ub- In addition, existing 
turbulent velocity and conformation tensor fields were 
restarted for computations where We c or L p was modi- 
fied. In these cases, the flow undergoes a transient time, 
where again sufficient statistics were collected after reach- 
ing a stationary state. 



According to Eqs. <j3j)-([5]) , the four dimensionless 
groups that can fully characterise the velocity and the 
conformation tensor fields are We c , L p , j3 and Re c , and 
they are tabulated below. ?he reasons behind the choice 
of the particular parameter values is outlined below. The 
rationale here follows the thorough parametric study by 
Li et al. jsj. 



TABLE I: Parameters for the DNS of viscoelastic turbulent channel flow. The friction Weissenberg number is defined by 



We T0 = 


T P U 2 T JV. 


LDR cases: 


A, B, D2, I, J; 


HDR cases: 


C, D, Dl, E, F, G, K; MDR case: H. 






Case 


We c 


We T0 




/3 


Re c 


Re-r 


L x x Ly x L z 


X N y x N z 


DR(%) 


Nl 








1 


2750 


123.8 


6.5tv5 x 25 x 1.5n6 


200 x 65 x 100 





N2 








1 


4250 


181 


4.5-7r<5 x 25 x tyS 


200 x 97 x 100 





N3 








1 


10400 


392.6 


2tt<5 x 25 x 0.5tt<5 


200 x 193 x 100 





A 


2 


15.4 


120 


0.9 


4250 


167.7 


4.5nS x 28 x nS 


200 x 97 x 100 


-14.2 


B 


4 


30.8 


120 


0.9 


4250 


147.3 


4.5n5 x 28 x n5 


200 x 97 x 100 


-33.8 


C 


7 


54 


120 


0.9 


4250 


121.8 


4.5n5 x 28 x n5 


200 x 97 x 100 


-54.7 


D 


9 


69.4 


120 


0.9 


4250 


118.3 


4.5n5 x 28 x n5 


200 x 97 x 100 


-57.3 


Dl 


9 


69.4 


60 


0.9 


4250 


124.7 


4.5nS x 25 x n5 


200 x 97 x 100 


-52.5 


D2 


9 


69.4 


30 


0.9 


4250 


150.3 


4.5n5 x 25 x n5 


200 x 97 x 100 


-31 


E 


11 


84.8 


120 


0.9 


4250 


113.3 


4.5n5 x 25 x ty5 


200 x 97 x 100 


-60.8 


F 


13 


100.2 


120 


0.9 


4250 


112.4 


4.57t<5 x 25 x tt5 


200 x 97 x 100 


-61.4 


G 


15 


115.6 


120 


0.9 


4250 


111.4 


4.57r<5 x 25 x ttS 


200 x 97 x 100 


-62.1 


H 


17 


131 


120 


0.9 


4250 


107.8 


8ttS x 25 x ttS 


200 x 97 x 100 


-64.5 


I 


2 


29.6 


120 


0.9 


10400 


323.3 


2-kS x 25 x 0.5tt5 


200 x 193 x 100 


-32.2 


J 


4 


22.3 


120 


0.9 


2750 


106.9 


6.5n5 x 25 x 1.5x6 


200 x 65 x 100 


-25.4 


K 


7 


39 


120 


0.9 


2750 


91.1 


6.5n6 x 25 x 1.571-5 


200 x 65 x 100 


-45.9 



Drag reduction effects are expected to be stronger at 
high Weissenberg numbers. In fact, higher levels of per- 
centage drag reduction at MDR have also been measured 
for higher Reynolds numbers 0, showing the Reynolds 
number dependence on drag reduction amplitude. There- 
fore, in this work, an extensive parametric study has been 
carried out by mainly varying We c for the computation- 
ally affordable Re c = 4250 to determine the impact of 
polymer dynamics on the extent of drag reduction. The 
Reynolds numbers considered here, Re c = 2750, 4250 and 
10400, are small in comparison to most experimental 
studies but fall within the range of most DNS studies of 
polymer induced turbulent drag reduction. Nevertheless, 
these Reynolds numbers arc sufficiently large for the flow 



to be always turbulent and allow to study the dynamics 
of viscoelastic turbulence. Different maximum dumbbell 
lengths were also taken into account to check their effects 
for the same We c and Re c . The chosen L p = b + 3 values 
are representative of real polymer molecule lengths which 
can be related through b « Nc/(vlfN) where Nc is the 
number of carbon atoms in the backbone of the polymer 
macromolcculc, cr s f is an empirical steric factor, N is the 
number of monomers and b is the finite dumbbell exten- 
sibility and is a large number [l4| . Note that in the limit 
b — y 00 , the Hookean spring-force law is recovered, which 
governs a linear spring. 

Low P values were used in most prior DNS to achieve 
high levels of drag reduction, in view of the attenuation of 
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the polymer-turbulence interactions due to the additional 
artificial diffusion term in the FENE-P model and their 
moderate Reynolds numbers, usually Re T = 5/8 v < 395. 
In fact, values as low as /3 = 0.4 have been applied thus 
amplifying viscoelastic effects so as to reach the HDR 
regime [32]. However, such low j3 values may lead to 
significant shear-thinning 1 [ll[ unlike in experiments of 
polymer drag reduction. The fact that the numerical 
scheme applied in our study for the FENE-P model is 
expected to capture the strong polymer-turbulence in- 
teractions allows the value of /3, which is inversely pro- 
portional to the polymer concentration, to be high, i.e. 
j3 = 0.9, representative of dilute polymer solutions used 
in experiments. 

The box sizes L x x L y x L z , where subscripts indi- 
cate the three Cartesian co-ordinates, were chosen with 
reference to the systematic study by Li et al. [3l[ of 
how the domain size influences the numerical accuracy. 
Specifically, they point out that long boxes are required 
in DNS of polymer drag reduction, particularly in the 
streamwisc direction because of longer streamwise corre- 
lations at higherpercentage DR, as opposed to the min- 
imal flow unit [33l | used in many earlier works. Differ- 
ent grid resolutions N x x N y x N z were tested for con- 
vergence. In particular, the following set of resolutions 
128 X 65 x 64, 200 x 97 x 100 and 256 x 129 X 128 were tried 
for Re To ~ 180 with the two latter giving identical mean 
velocity profiles and not significantly different rms veloc- 
ity and vorticity profiles. Similarly, grid sensitivity tests 
were carried out for the other Re To cases. Eventually, the 
resolutions for each Newtonian fluid computation were 
validated against previous ly p ublished databases for the 
corresponding Re ro cases [34M36j |. Note that if the res- 
olutions for Newtonian turbulent computations are ad- 
equately resolving the flow scales, then the same reso- 
lutions arc sufficient for viscoelastic turbulent computa- 
tions, since the size of vortex filaments in these flows 
increases while their number decreases as drag reduces 

For a given resolution, viscoelastic computations re- 
quire approximately 4 times more memory and 2 times 
more CPU time per time step compared to the Newto- 
nian case. The time step At used in viscoelastic compu- 
tations is typically a factor of 5 smaller than that used 
in the Newtonian cases due to the stricter CFL condition 
of the present numerical method for the FENE-P model 
(see Eq. (|A.9|) in the appendix and [26[ for more details 
on the time step constraint using compact schemes). Ul- 
timately, the viscoelastic computations require approxi- 
mately 10 times more CPU resources than the Newtonian 
computations for a given computational time period. 



III. VISCOELASTIC TURBULENCE 
STATISTICS 

A. Polymer drag reduction 

Since the computations are performed with a con- 
stant flow rate by adjusting the pressure gradient, DR 
is manifested via a decrease in skin friction, i.e. lower 
Rc r = 5/6 v values as drag reduces. Here, we define per- 
centage drag reduction as a negative quantity 

d 0>> _ f^^£l\ I 2 21 

dx [ dx I „ K -U%\ 

DR = -±- s J — x 100% = n , 10 x 100% 

= ((f^) 2 - 1 )*™ % (12) 

with u 2 . = — - . Variables with and without subscript 
in Eq. (fT2"j) refer to Newtonian 2 and viscoelastic fluid 
flow, respectively. Note that for a direct comparison be- 
tween the various cases with different skin friction we 
choose our plots to be presented in terms of y/S instead 
of y/$v because 6 is the same for all our computations, 
whereas 8 V increases with drag reduction. 

Figure [1] depicts the capability of the current numer- 
ical scheme for the FENE-P model to enable stronger 
polymer-turbulence interactions than artificial diffusion 
methods. Higher values of percentage drag reduction as 
function of Weissenberg number are obtained comparing 
with earlier DNS studies (see Fig. lb in (37} ) without 
the need for low /3 values [32[. These DR values extend 
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FIG. 1: (Colour online) Variation of percentage drag reduc- 
tion with Weissenberg number. 



2 Any departure from the Newtonian behaviour, i.e. <r;j oc Sij, 
with some constant of proportionality independent of the rate of 
the shear stress increases slower than linear 012 oc S12 strain, could be called non-Newtonian. 
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throughout the drag reduction regimes. The MDR limit 
is approached in this case at |DR| ~ 65% because of 
the moderate Re c in our computations. Even so, this 
amount of drag reduction falls within the MDR regime, 
based on the classification of drag reduction by Warholic 
et al. Q, allowing to study the MDR dynamics of the 
polymer molecules and their effects on the flow in this 
asymptotic state. 



B. Effects of polymer extensibility and Reynolds 
number 
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The effects of maximum dumbbell extensibility is 
briefly considered for three different extensibilities L p = 
30, 60 and 120 but the same We c and Re c (see D cases 
in Table HJ . Figure [1] shows that the extent of drag re- 
duction is amplified by longer polymer chains consistent 
with other DNS studies J3J, [38( . This effect is related 
to the fact that the average actual length of the dumb- 
bells, represented by the trace of the conformation tensor 
(Cfcfc), increases further for larger L p according to Fig. 
[2^,, inducing stronger influence of the polymers on the 
flow. The percentage increase, however, of the polymers 
extension is less for larger FENE-P dumbbells (see Fig. 
[2Jd) , suggesting that large polymer molecules could be less 
susceptible to chain scission degradation, which causes 
loss of drag reduction in experiments 0] . The near- wall 
turbulence dynamics play an important role for all three 
cases, as most of the stretching happens near the wall, 
where the highest fluctuating strain rates are expected. 
Eventually, the largest maximum length, i.e. L p = 120, 
was used for the rest of the computations considered in 
this work in order to explore the polymer dynamics at 
effective drag reductions, which are interesting not only 
fundamentally but also in many real life applications. 

Based on DNS with artificial diffusion methodology, 
Housiadas and Beris [39[ claim that the extent of drag 
reduction is rather insensitive to Reynolds numbers in 
the range between 125 < Re To < 590 for LDR flows. 
On the other hand, avoiding the use of artificial diffusion 
in our study, the Reynolds number dependence on drag 
reduction for cases with identical We c values but different 
Reynolds numbers, i.e. Re c = 2750,4250 and 10400, is 
obvious by comparing DR of case A with I and case B 
with J (LDR regime), as well as case C with K (HDR 
regime), where the percentage DR increases for higher 
Re c at all instances (see Tabic [I]). This Reynolds number 
dependence is further depicted in the polymer dynamics 
through the profiles of (Ckk)/L p , which amplify closer to 
the wall due to more intense strain rates in this region at 
higher Re c and collapse towards the centre of the channel 
(see Fig. [3]). The disparate behaviour of (Ckk) with 
respect to y/8 due to the Reynolds number dependence 
is anticipated by the broader spectra of flow time scales 
that are encountered at higher Re c by the dumbbells with 
fixed relaxation time scale. The fact that the current 
DNS could capture the Reynolds number dependence on 
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FIG. 2: (Colour online) Effect of maximum dumbbell length. 
Plots of (a) average actual dumbbell extensibility (Ckk) and 
(b) percentage average dumbbell extensibility (Ckk)/L p as 
functions of y/8. Note: case D (L p = 120); case Dl (L p = 60); 
case D2 (L p = 30). 



drag reduction and polymer dynamics, emphasises once 
more the strong polymer-turbulence interactions that can 
be captured by the present numerical approach even at 
low levels of drag reduction. 

It is essential to note at this point that the interme- 
diate dynamics between the von Karman and the MDR 
law, i.e. the LDR and HDR regimes, are non-universal 
because they depend on polymer concentration, chemical 
characteristics of polymers, Reynolds number, etc. 
Here, this is illustrated by the maximum dumbbell length 
and Reynolds number dependencies of the polymer dy- 
namics in Figs. [2] and [3l respectively. However, at the 
MDR limit, which is achieved at We c ;§> 1 and Re c ^> 1, 
the dynamics are known to be universal 0, H[, i.e. inde- 
pendent of polymer and flow conditions. 
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Case B 
Case C 
Case I 
Case J 
Case K 





FIG. 3: (Colour online) Effect of Reynolds number on per- 
centage average dumbbell extensibility as function of y/5. 
Identical symbols correspond to cases with the same We c val- 
ues. Note: Compare case A (We c = 2, Re c = 4250) with case 
I (We c = 2, Re c = 10400); case B (We c = 4, Re c = 4250) with 
case J (We c = 4, Re c = 2750); case C (We c = 7, Re c = 4250) 
with case K (We c = 7, Re c = 2750). 



FIG. 4: (Colour online) Mean velocity profiles versus y+ for 

the LDR, HDR and MDR regimes. : U+ = y+, - - -: U+ = 

^ \ogy+ + 6.0, • • •: U+ = j±j logy + - 17. Note: case N2 
(DR = 0%); case A (DR = -14.2%); case B (DR = -33.8%); 
case D (DR = -57.3%); case G (DR = -62.1%); case H 
(DR = —64.5%). 



C. Mean and fluctuating velocity statistics 

The picture of drag reduction can be analysed in fur- 
ther detail with the statistics of the turbulent velocity 
field introduced in Figs. 0] and [5j The distinct differ- 
ences in the statistical trends of the turbulent velocity 
field between the LDR and HDR regime, that have been 
observed experimentally, arc clearly identified in these 
results. For clarity, a few indicative cases from the data 
of TableUhave been chosen for plotting, representing the 
LDR, HDR and MDR regimes for different Weissenberg 
numbers at Re c = 4250. 

According to Fig. @] and noting that j3 = 0.9 for all 
viscoelastic cases, all mean velocity profiles collapse in 
the viscous sublayer y + < 10 to the linear variation U+ = 
f3~ 1 y +1 which can be deduced for viscoelastic flows, by 
rewriting Eq. (jlip in viscous scales 



dU+ (uV) | (<r 12 ) 
dy+ u% v? r 



and neglecting the normalised Reynolds and mean poly- 
mer shear stress in the viscous sublayer y + — > (see 
also section |V| . Figure |4] presents the clear impact of 
percentage DR on the mean flow with the skin friction 
decreasing and the mean velocity increasing away from 
the wall in comparison to the Newtonian case N2 as a 
result of higher We c values at the same Re c . The profile 
of the Newtonian case N2 is in agreement with the von 
Karman law Eq. (JTJ, which does not hold for viscoelastic 
turbulent flows. Specifically, the curves of cases A and 
B (LDR regime) are shifted upwards with higher values 
of the intercept constant B, i.e. parallel to the profile 



of the Newtonian flow (see Fig. [4j, increasing DR. This 
picture is consistent with the phenomenological descrip- 
tion by Lumley d, |40| . where the upward shift of the 
inertial sublayer can be interpreted MS cl thickening of 
the buffer or elastic layer for viscoelastic flows, which is 
equivalent to drag reduction. HDR cases D and G exhibit 
different statistical behaviour than LDR flows with the 
slope of the log-region increasing until the MDR asymp- 
tote is reached by case H. Overall, the same behaviour 
across the extent of drag reduction in viscoelastic tur- 
bulent flows have been seen in several experimental and 
numerical results @. 

Different statistical trends between low and high d rag 
reduction have also been observed experimentally [fl 41] 
for the rms streamwise velocity fluctuations normalised 
with u T . Figure [5^ illustrates the growth of the peak 
in the profile of u' + for LDR case A and B at low 
We c and a notable decrease for the rest of the cases 
at HDR/MDR with high Wc c values. The peaks move 
monotonically away from the wall throughout the drag 
reduction regimes indicating the thickening of the elas- 
tic layer, which is compatible with the behaviour of the 
mean velocity profile. 

Note that this is the first time that a DNS computation 
can so distinctly attain this behaviour. This is attributed 
to the accurate shock-capturing numerical scheme we ap- 
plied for the FENE-P model in this study. It has to 
be mentioned however that there have been three earlier 
studies [HI, |37], IHJ , which use the artificial diffusion algo- 
rithms for FENE-P and showed similar but not as clear 
trends for u'< in a DNS of viscoelastic turbulent channel 
flow. In fact, Min et al. j37| reached the HDR/MDR 
regime at roughly |DR| ~ 40%, clearly very low to afford 
the correct dynamics and Ptasinski et al. [32j had to 
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FIG. 5: (Colour online) Rms velocity components for the 
LDR, HDR and MDR regimes, (a) Streamwise u' + , (b) wall- 
normal v' + and (c) spanwise w' + profiles versus y/S. Note: 
case N2 (DR = 0%); case A (DR = -14.2%); case B (DR = 
-33.8%); case D (DR = -57.3%); case G (DR = -62.1%); 
case H (DR = -64.5%). 



use j3 = 0.4 to approach HDR/MDR, encountering con- 
siderable shear-thinning effects. It is interesting to men- 



tion that other recent studies [3l], |43| , using the artificial 
diffusion methodology, with more extensive Weissenberg 
number data and high /3 values, have not been able to ob- 
tain this transition effect on the statistics of u\ between 
the drag reduction regimes. 

Finally, the wall-normal v' + and spanwise w', rms ve- 
locity fluctuations in Figs. [Bp and [Bp, respectively, are 
continuously attenuated while DR is enhanced by in- 
creasing the polymer relaxation time scale. Again, the 
monotonic displacement of their peaks towards the cen- 
tre of the channel as drag reduction amplifies is consistent 
with that of the mean velocity profile and with experi- 
mental and other numerical studies 21. 



D. Fluctuating vorticity statistics 

The rms statistics of the fluctuating vorticity field nor- 
malised by viscous scales, i.e. u>' + = oj'S v /u t , are pre- 
sented in Fig. [Bj for representative cases from Table U at 
various levels of drag reduction. The streamwise vortic- 
ity fluctuations u' x demonstrate a persistent attenuation 
along the normalised distance y/S as drag reduction en- 
hances due to the increase of We c (see Fig. [BJi). In the 
near-wall region y/S < 0.2 of Fig. [BJi there is a character- 
istic local minimum and maximum that could be inter- 
preted as corresponding to the average ed ge a nd centre of 
the streamwise vortices, respectively |3ll |44| . Then, the 
average size of these large streamwise vortices is roughly 
equal to the distance between these two peaks. The fact 
that these peaks are displaced away from each other and 
at the same time away from the wall, as DR builds up, 
implies an increase in the average size of the streamwise 
vortices and a thickening of the buffer layer, respectively, 
in agreement with earlier works [U, HH, Hfl Hg . The atten- 
uation in the intensity of uj' x+ provides evidence for a drag 
reduction mechanism based on the suppression of the 
near- wall counter-rotating steamwise vortices [46|, l47j . 
which underpin considerable amount of the turbulence 
production [48j |. 

The wall-normal rms vorticity is zero at the wall due 
to the no-slip boundary condition and reaches its peak 
within the buffer layer (see Fig. [Bp). The intensity of 
uj' y+ is reduced for all levels of drag reduction according 
to Fig. [6p, with the position of the near-wall peaks mov- 
ing towards the centre of the channel as We c becomes 
larger, representing once more the thickening of the elas- 
tic layer in a consistent way. Most of the inhibition of 
u>' happens near the wall and slightly towards the cen- 
tre of the channel only for the HDR/MDR cases G and 
H, i.e. for |DR| > 60%. 

Figure [Bp shows a more interesting behaviour for u>' z , 
where the spanwise vorticity fluctuations decrease in 
the near- wall region y/S < 0.2 and increase further 
away while drag reduces. This effect may be related 
to the transitional behaviour of u' + between the LDR 
and HDR/MDR regimes (see Fig. [B^) plus the contin- 
uous drop of v' + (see Fig. [Bp) in viscoclastic drag re- 
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FIG. 6: (Colour online) Rms vorticity profiles for the LDR, 
HDR and MDR regimes, (a) Streamwise u)' x+ , (b) wall-normal 
oj' y+ and (c) spanwise u' z profiles versus y/8. Note: case N2 
(DR = 0%); case A (DR = -14.2%); case B (DR = -33.8%); 
case D (DR = -57.3%); case G (DR = -62.1%); case H 
(DR = —64.5%). 

duced flows. As a final note, ui' > uL > u/„ in 

■**+ y+ 

the viscous sublayer, i.e. y/8 < 0.05 for all cases and 



uj' z+ ~ oj' x+ ~ Lo' in the inertial and outer layer for the 
Newtonian case N2. However, uj' z+ > uj' y+ > uj' x+ away 
from the wall when drag reduces for viscoelastic flows, 
which manifests the dominance of small scale anisotropy 
in the inertial and outer layer at HDR and MDR. 

IV. CONFORMATION AND POLYMER STRESS 
TENSOR 

Before looking at the mean momentum and energy 
balance, the study of the conformation tensor field is 
essential to get an understanding of the polymer dy- 
namics in support of the results provided by this new 
numerical method for the FENE-P model in turbulent 
channel flow. The symmetries in the flow geometry de- 
termine properties of tensor components in the average 
sense [49|. In the current DNS of turbulent channel 
flow, statistics are independent of the z direction and 
the flow is also statistically invariant under reflections 
of the z co-ordinate axis. Therefore, for the probability 
density function f(Q;x,t) of a vector Q, these two con- 
ditions imply df/dz = and f(Qi,Q2,Qz',x,y,z,t) = 
f(QiiQ2,—Q3',x,y,—z,t). Then, at z = reflectional 
symmetry suggests that (Q3) = — (Q3) => (Q3) = and 
similarly for (Q1Q3) = (Q2Q3) = 0. So, in this case the 
mean conformation tensor reduces to 

f(Cu) (Cia) \ 
(Cij) = (C 12 ) (C22) (14) 
V (C 33 )J 

where the non-zero components scaled with L p are 
presented in Figs. [7] and [8] with respect to y/S 
for cases at various drag reduction regimes (see Ta- 
ble |T|. The zero components in our study have 
been found to be zero within the precision accuracy. 
Turbulent channel flow is also statistically symmet- 
ric about the plane y = 8. Therefore, this re- 
flectional symmetry imposes f{Q\,Q2,Q3]x,y,z,t) = 
f{QiT~Q2,Q3\x,—y,z,t), which implies that the nor- 
mal components of (Cij) are even functions and (C12) is 
an odd function comparable to the Reynolds stress tensor 
components. 

The normalised trace of the mean conformation tensor 
(C k k)/Lp is plotted in Fig. [T^i together with {Cn)/L 2 p . 
Notice that the dominant contribution in the trace comes 
from (Cn), i.e. (Ckk) — (Cn) at all Weissenberg num- 
bers, reflecting on average a strong preferential orien- 
tation of the stretched dumbbells along the streamwise 
direction. The fact that {C n ) > (C 12 ) ^ (C* 33 > > (C 22 ) 
denotes the strong anisotropic behaviour of the mean 
conformation tensor caused by the mean shear in turbu- 
lent channel flow. This anisotropy influences the statis- 
tics of the fluctuating velocity field particularly at small 
scales, as was mentioned in section UlIDI The curves of 
(Cn)/Lp and consequently of (Ckk)/L p constantly rise 
with most of the stretching happening close to the wall 
and growing towards the centre of the channel, since 
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FIG. 7: (Colour online) Profiles of (a) {C kk )/L 2 p (line- 
symbols) and (Cn)/Lp (solid lines) , (b) (Ci2)/Lp as functions 
of y/8 for the LDR, HDR and MDR regimes. Note: case A 
(DR = -14.2%); case B (DR = -33.8%); case D (DR = 
-57.3%); case G (DR = -62.1%); case H (DR = -64.5%). 



higher values of polymer time scale are influenced from a 
wider spectrum of flow time scales. A local minimum and 
a maximum emerge in the near- wall re gion y/8 < 0.2, in- 
duced by the streamwise vortices (42l. |50|. These peaks 
move apart from each other and away from the wall for 
higher We c values. Figure [7^, also shows that the ampli- 
tudes of these peaks seem inversely proportional to the 
peak amplitudes of oj' x+ as drag reduces (see also Fig. 
Eli). 

Moreover, as We c increases the profiles of {C\2)/L 2 
and (C33)/Lp amplify, reaching their peaks at not much 
different y/6 for each We c case (see Figs. [7p and[Sp). 
In particular, the values of {C\2)/L 2 at the wall are de- 
pendent on the polymer relaxation time scale unlike for 
{Cz?,)/L 2 p . On the other hand, the values of (C 33 )/Lp 
depend on Weisscnberg number at y = 5 in contrast to 
(Ci2)/Lp, which is zero for all cases because of the sym- 
metry mentioned earlier. The behaviour of (C22)/Lp in 
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FIG. 8: (Colour online) Profiles of (a) (C 22 ) / 'L% and (b) 
<C 3 3>/ip as functions of y/8 for the LDR, HDR and MDR 
regimes. Note: case A (DR = -14.2%); case B (DR = 
-33.8%); case D (DR = -57.3%); case G (DR = -62.1%); 
case H (DR = -64.5%). 

Fig. [8^ is more peculiar with respect to We c , with the 
profiles increasing within the LDR regime and attenuate 
for HDR and MDR cases, in a similar manner to u' + (see 
Fig. [5^). Its peak values are achieved closer to the core 
of the channel in comparison to the rest of the confor- 
mation tensor components. This points out the different 
flow time scales that are important for (C22), exemplify- 
ing the complex dynamics of the polymers, even in this 
simple mechanical model. 

It is interesting to mention that the components of 
(Cij) have different asymptotic rates of convergence to- 
wards the limit We c —toe. It is known that for We c 3> 1 
the upper bound for the trace is (Ckk) < L 2 and subse- 
quently in our case (Cn) < L 2 (see Fig. [7^), where this 
upper bound is far from achieved in our computations. 
This result demonstrates that highly stretched polymers 
are not required for the manifestation of drag reduc- 
tion or even of the MDR asymptote, as de Gennes [5l[ 
claims against Lumley's [|| assumption of a coil-stretch 
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transition, i.e. highly stretched polymer molecules, for 
the enhancement of intrinsic viscosity. The components 
(C12 ) /Lp and (C33 ) / l? v seem to have almost reached their 
asymptotic limit with the MDR case H according to Figs. 
[7p and [8p, respectively. Finally, (C22) / Lp has not yet 
converged to its limit, decreasing with a slow rate to- 
wards very small values for high We c . In fact, it has 
been argued theoretically that (C22) — > in the limit of 
infinite Weissenberg number 0, |52| | . 

Polymer stresses are non-linear with respect to the con- 
formation tensor and their asymptotic scaling with Weis- 
senberg number is a key element for the understanding of 
the polymer dynamics at MDR. Hence, following Benzi 



et at [53[ consider the FENE-P model integrated over 



the x, z spatial directions and time, assuming statistical 
stationarity and homogeneity in x and z 

(u2d X2 Cij) = (C ik d Xk Uj) + (C jk d Xk Ui) 

-^-(fiCkkPy-Sij). (15) 

Then, taking the Reynolds decomposition of the velocity 
field Ui = (ui) + u[, one obtains 

1 



We, 



■{Steven - Sij) = (C^dM+iC^d^+Qij 

(16) 

where = (CiA^) + (C^c^u-) - (u 2 d X2 C tJ ). 

Therefore, the average polymer stress tensor defined by 
Eqs. (01 takes the form 



1-/3 
Rc c 



Qn 
Q12 



Q22 



Q 



12 



Q 



13 



Q 



2:; 




Now, the important assumption at the limit of a local 
Weissenberg number Wes = t p j-(u) — > 00 is that Qn 
and Q12 can be neglected, considering the polymers to 
be stiff, i.e. (7y — > (C,j), mostly in the main stretching 
directions and the correlations between fluctuating con- 
formation tensor and velocity fields in the other Cartesian 
directions to remain minimal at this limit. In this case, 
as a result 

(a u ) = A 1 ^-2(C 12 )d X2 <« x ) (18) 
Re, 



{ai2)=A 2 1 —^{C22)d X2 {u l ) 
Re, 



(19) 



where A\ and A2 are expected to be independent of y 
and equal to 1 at some intermediate region in the flow as 
Wes 3> 1- This hypothesis is checked in Fig. [S] against 
various viscoelastic DNS from Table Q] 

Figure [SJi shows clearly that A\ tends to a constant 
and reaches 1 in the region < y/S < 0.8 for high Wes 
values justifying that Qn can be neglected for HDR and 




FIG. 9: (Colour online) Scalings of the compensated polymer 
stress components (a) A x = (an)/ (-g^f 2(Ci 2 }^-) and (b) 

A2 = (C12)/ (^^r (C22) with respect to y/S. Note: case 

A (DR = -14.2%); case B (DR = -33.8%); case D (DR = 
-57.3%); case G (DR = -62.1%); case H (DR = -64.5%). 



MDR cases. Note that A\ deviates from 1 towards the 
centre of the channel because Wes becomes small in this 
region. A2 is approximately independent of y in some 
intermediate region for almost all cases and appears to 
tend towards 1 as Wes increases (see Fig. [9p). However, 
the polymer relaxation time scales used in this study are 
not sufficiently high for A<2 — > 1. So, for our DNS re- 
sults the polymer shear stress can be considered to be 
M (x ^{C 22 )d X2 { Ul ) in a range 0.2 < y/S < 0.7. 
It is appealing to see that (C22) is the component in- 
volved in the MDR dynamics, bearing in mind that 
(C u ) > (Cia) ^ (C 33 ) > (C22). In the end, both Figs. 
[TJj and [9^i confirm the claims that (C12) has reached its 
asymptotic limit within the Weissenberg numbers con- 
sidered at this particular Reynolds number in this study, 
unlike (C 22 ) (see Figs. UK and Hp). 
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V. SHEAR STRESS BALANCE 

The balance of total shear stress Eq. (fTTj) is con- 
sidered in this section. The total shear stress in vis- 
coelastic turbulent channel flow contains the shear stress 
(3v-^(u) coming from the mean flow, the Reynolds stress 

— {u'v') rising from turbulence and the mean polymer 
shear stress (C12) due to polymers in the flow, which 
is also referred to as the Reynolds stress deficit since 
v ~ckj( u ) ~ (u'v 1 ) 7^ v% (1 — y/5). In addition, the gradient 
of the Reynolds shear stress rises due to the non-linear 
term in the Navier-Stokes equations and it is related to 
the lift force (Magnus effect) experienced by a vortex line 
exposed to a velocity u [5J]. The Reynolds shear stress 
and the rotational form u x us of the non-linear term 
in Eq. ([3]) are related through the following equation 
J^(u'v') = (v'oj' z ) — (w'uj' y ). It is true that tangles of 
very intense and slender vortex filaments can exist down 
to very fine scales creating an energy drain on the mean 
flow. These tangles can be seen as one form of inter- 
mittency in turbulent flows (55j . However, it is unclear 
how the intermittency of the vorticify field u> affects the 
averages (v'ui' z ) and (w' ' ui' y ) and thereby Reynolds shear 
stress via the relation -j^(u'v') = {v'lo' z ) — (w'ui' y ). 

The viscous stress of the solvent, the Reynolds shear 
stress and the mean polymer shear stress normalised with 
viscous scales are presented in Fig. [TU] at different levels 
of percentage DR for cases with the same Re c from Table 
HI At the wall, the no-slip boundary condition enforces 

— (u'v')\y_ = 0. The wall shear stress is governed by 

90% viscous as well as 10% polymer contribution for all 
viscoelastic cases as opposed to the Newtonian case N2. 
Viscosity is the dominant parameter in the near- wall re- 
gion but becomes more influential in the outer regions 
as drag reduction enhances. This is clear from Fig. 110b , 
where /3-^-U+ increases monotonically towards the cen- 
tre of the channel as Wc c increases. Viscoelastic effects 
become also more significant towards the centre of the 
channel for higher We c cases. Reynolds shear stress, on 
the other hand, is constantly decorrelated at higher per- 
centage DR with its peak shifting away from the wall. 
It is interesting to see that for the HDR case D —(u'v') 
and ((T12) are comparable and as MDR is approached the 
polymer shear stress plays an increasingly fundamental 
role in sustaining turbulence due to the vast attenua- 
tion of the Reynolds shear stress at these finite Reynolds 
number computations. This becomes apparent in the 
next section by analysing the turbulent kinetic energy 
budget. 

Notice that Reynolds shear stress remains finite at 
MDR confirming the experimental measurements by 
Ptasinski et al. [4l[ against the complete depletion of 

— {u'v') reported by Warholic et al. 0] and their sub- 
sequent claim that turbulence is sustained entirely by 
polymer stresses. What can be said theoretically on this 
controversy is the following. Consider first the limit of 
Wes — > 00, where A2 — > 1 for Eq. (jTSJ) even at finite 
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FIG. 10: (Colour online) Profiles of (a) viscous shear stress, 
(b) Reynolds shear stress and (c) mean polymer shear stress 
versus y/5 for the LDR, HDR and MDR regimes. Note: case 
N2 (DR = 0%); case A (DR = -14.2%); case B (DR = 
-33.8%); case D (DR = -57.3%); case G (DR = -62.1%); 
case H (DR = -64.5%). 



Reynolds numbers, as Fig. [!j> suggested. Then, the to- 
tal shear stress balance Eq. (fTTj) can be rewritten using 
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Eq. CEU) 



v(P + (1 - P)(C 22 )) 



dy 



(«V) 



I) (20) 



where v e ff(y) = v{P + (1 — /3)(C22)) is an effective vis- 
cosity similar to the one encountered in Lumley's phe- 
nomenology [H, Eo[. Now, when Wes 3> 1 assume that 
(C22) becomes minimal based on theoretical claims by 
|9|, |52| and observational indications in this study. Then, 
for high enough Reynolds number along the universal 
MDR asymptotic line, i.e. taking first the infinite Weis- 
senberg number limit and then the infinite Reynolds 
number limit, one might expect an intermediate region 
5 V <C y -C 5 of approximately constant Reynolds shear 
stress, i.e. —(u'v')/u^. — > 1, implied by Eq. ([2T)|) when 
taking the limits of y/5 — > and y/S v — > 00 with the rea- 
sonable assumption that vj3-^{u) — > as y 5 U . This 
statement suggests that the classical way of turbulence 
production does not vanish in the infinite Weissenberg 
and Reynolds number limit. 

Ultimately, the conjecture here is that (012) becomes 
negligible in the stress balance Eq. (1TTT) when carefully 
taking the double Weissenberg and Reynolds number 
limit in the right order, so that we are on the universal 
MDR asymptotic line. This, however, does not indicate 
that drag reduction is depleted, it rather suggests that 
the MDR asymptote could be entirely determined by the 
energetics in these infinite limits. Nevertheless, polymers 
play a crucial role in the dynamics at MDR and this will 
be explored further in the next section. 



VI. POLYMER-TURBULENCE DYNAMICAL 
INTERACTIONS 

The balance equation for the turbulent kinetic energy 
of a viscoelastic fluid provides further insight into the dy- 
namical interactions between polymers and turbulence. 
Assuming statistical stationarity and homogeneity in x 
and z directions for the mean turbulent kinetic energy 
balance and integrating over the y direction, we obtain 



Vdy = e N dy+ £pdy 



(21) 



with no contribution from the transport terms due to the 
no-slip boundary condition, using the divergence theo- 
rem. The turbulence production by Reynolds shear stress 
is denoted here by V = -(uV)^(u), the viscous dissi- 
pation rate En = 2^/3 (sijSij) and the viscoelastic dissipa- 
tion rate Ep = {c'ijdxjU'i), which arises due to fluctuating 
polymer stresses. Note that Ep has a dual nature, i.e. it 
can serve either as dissipation or production depending 
on the signs of the polymer stress fluctuations and that 
of the fluctuating velocity gradients. 

Figure QT] presents each term of Eq. (|2"Tj) normalised 
by 5 v /u^ with respect to We T0 for all cases from Ta- 
bic U at Re c = 4250. An asymptotic behaviour to a 
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FIG. 11: (Colour online) Terms of the y-integrated turbulent 
energy balance with respect to We T0 . 



marginal flow state can be observed by increasing the 
polymer relaxation time scale with a vast attenuation oc- 
curring in the total production and viscous dissipation, 
while viscoelastic dissipation grows mildly in the LDR 
regime and constantly decays within HDR and MDR. 
Overall, J Ep dy becomes pivotal in the dynamics of the 
flow relative to J V dy and J En dy for HDR and MDR 
flows. Most importantly J Epdy < for high We ro val- 
ues according to Fig. 1111 in agreement with experimen- 
tal measurements [4l|, implying that polymers somehow 
can sustain turbulence by producing turbulent kinetic 
energy. Notice, that in this plot both dissipations are 
presented as positive quantities and this was done on 
purpose to emphasise the interplay between production 
and viscous dissipation from LDR to HDR. It is notewor- 
thy that J V dy > J En dy for LDR cases A and B but 
j V dy < J En dy for HDR cases and gets even smaller 
as drag reduction approaches its maximum limit. This 
observation hints that polymer dynamics get somehow 
involved in the production of turbulent kinetic energy so 
that turbulence does not die out at HDR and MDR. 

Lets now look in more detail at the profiles oiV, £n 
and Ep scaled by 8 v /v%. with respect to normalised dis- 
tance from the wall y/5 for representative cases at various 
levels of drag reduction from Table|T](see Fig. [T^]). Dissi- 
pation represents drain of energy, hence, En and Ep have 
been plotted here as negative quantities. 

The production of turbulent energy by Reynolds 
stresses, which is continuously reduced over the extend of 
drag reduction as a function of We c , serves to exchange 
kinetic energy between the mean flow and the turbulence. 
The local peak of V is reached within the buffer layer 
and in fact for Newtonian flows we can easily show that 
the maximum production occurs where —(u'v') = f-^{u) 

and Vmax5 u /u^. < | [4^ |. The peak turbulence produc- 
tion within the LDR regime also occurs at the intersec- 
tion point of viscous and Reynolds shear stress (compare 
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case H (DR = -64.5%). 



Figs. [TUk and [TUb with Fig. [T2"k). which shifts away 
from the wall as We c increases, indicating the thickening 



of the elastic layer. However, for HDR and MDR cases 
'PmaxSu/ U T is within 0.1 < y/8 < 0.3, where the maxi- 
mum Reynolds stress roughly appears, without following 
the —(u'v') = f3v-^{u) intersection point, which does not 
even exist for cases G and H (see Figs. [TUk and 1 10b). 

Viscous dissipation exhibits monotonic attenuation as 
drag reduces for higher values of We c with the maximum 
dissipation arising at the wall for the Newtonian case N2 
and the LDR cases A and B (see Fig. [12b). Although 
the kinetic energy is zero at the wall since u'\ y= o = im- 
posed by the no-slip boundary conditions, the fluctuating 
strain rate and consequently en is non-zero. At high per- 
centage DR, we surprisingly observe that the highest fluc- 
tuating strain rates are encountered away from the wall 
providing a completely different picture of the near-wall 
dissipation dynamics. The local kink in the buffer/clastic 
layer, which arises due to intense activity in this region, 
exists at corresponding y/8 with V max 8 v /u^ r for all cases 
considered in Fig. 112b and becomes a global minimum 
for the HDR and MDR cases, dominating the profiles of 
viscous dissipation. 

The profiles of viscoelastic dissipation obey a charac- 
teristic transitional trend similar to what has been al- 
ready observed for u' + (see Fig. [5k) and (C22) (see Fig. 
[Sk) from LDR to HDR regime, as We c increases. In 
detail, the curves of LDR cases A and B shift down- 
wards increasing viscoelastic dissipation but those of the 
HDR/MDR cases move upwards enhancing the positive 
nature of —ep8 v /v? T . The dual nature of ep is clearly de- 
picted in Fig. 112b with polymers dissipating and produc- 
ing turbulent kinetic energy in different regions, which 
depend on the polymer relaxation time scale at a given 
Reynolds number. A Reynolds number dependence of 
these regions is expected owing to the effect of different 
flow time scales on dumbbells with a particular relaxation 
time scale. Figure [T3] compares cases of identical Weis- 
senberg numbers and different Reynolds numbers (see 
Table [TJ) , illustrating a weaker Re c dependence on vis- 
coelastic dissipation in comparison to the stronger We c 
dependence in Fig. IT2"b . particularly at HDR and MDR. 
The part of the total dissipation that occurs in the three 
regions defined by the profile of viscoelastic dissipation in 
Fig. [T2"b can be estimated based on the profiles in Figs. 
RT2"b and [JJk. Approximately 15% to 25% of the total 
dissipation takes place in the first region, 25% to 60% in 
the second region and 20% to 70% in the third region. In 
other words, the majority of the total dissipation occurs 
away from the wall. 

Now, considering each component of the correlation 
matrix ep = {^ijd Xj u'/) , where summation applies over 
the indices i and j, we can observe that components with 
i = 2, 3 can be ignored, with most of the contribution as- 
cribed to i = 1 components according to Fig. [T31 which 
is very similar to Fig. [T2"b . The qualitative features 
of ep are clearly captured by (ay^uj), simplifying the 
underpining dynamics of viscoelastic dissipation. How- 
ever, to be precise, ep is neither exactly approximate 
nor proportional to (cr^ ■ d Xj v! x ) . Note that the positive 
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FIG. 13: (Colour online) Effect of Reynolds number on vis- 
coelastic dissipation as function of y/8. Identical symbols 
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FIG. 14: (Colour online) — (a' 1 jd x .u'i)8 v /u^. as function of y/8 
for the LDR, HDR and MDR regimes. Note: case A (DR = 
-14.2%); case B (DR = -33.8%); case D (DR = -57.3%); 
case G (DR = -62.1%); case H (DR = -64.5%). 



nature of ep is caused by the correlations — {a'^dx^u'^) 
and —(&i3d X3 u'i) (see Figs. [T5k andlTSb). The rest of the 
components are negative for all cases considered here and 
decrease monotonically as Wc c increases like — (&i2dx 2 u 'i) 
in Fig. 115b . The only exception though is — (c^c^itg), 
which also exhibits a dual trend, negligible however in 
comparison to the components presented in Fig. 1151 
Finally, the correlations in Figs. [l~5k and [l~5b are also 
responsible for the transitional behaviour of viscoelastic 
dissipation profiles from LDR to HDR discussed earlier. 
The current picture of the dual nature of ep was first 
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FIG. 15: (Colour online) Profiles of viscoelastic dis- 
sipation components for the LDR, HDR and MDR 
regimes. (a) — {a'^d x u')8 v /u^, (b) — (p' V xdy^)8 v /l£ r and 
(c) —(a' 13 d z u')8v/v% as function of y/8. Note: case N2 
(DR = 0%); case A (DR = -14.2%); case B (DR = -33.8%); 
case D (DR = -57.3%); case G (DR = -62.1%); case H 
(DR = -64.5%). 



predicted by Min et al. [56J at low Weissenberg num- 
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bers, adding an artificial diffusion term to numerically 
solve the FENE-P model. However, the present DNS 
are the first to capture so clearly these regions through- 
out the drag reduction regimes, predicting the appropri- 
ate dynamics at corresponding We c values. Once more, 
this is attributed to our numerical approach applied here 
for the FENE-P model that is able to capture stronger 
polymer-turbulence interactions than algorithms based 
on artificial diffusion. There are even results using the 
artificial diffusion methodologies that erroneously predict 
polymers never feeding energy back to the flow [9|, |32| . 
Hence, in view of the current distinctly transparent ob- 
servations a conceptual model for the mechanism of drag 
reduction is deduced in the next section. 



VII. DRAG REDUCTION MECHANISM 

The recent review on polymer drag reduction by White 
and Munghal (H reports that the numerical evidence is 
somewhat conflicting regarding the flow regions where 
polymers extend and contract. In this study, these re- 
gions can be identified by applying the Reynolds decom- 
positions Ui — (ui) + u\ and = (ov,- )+a J ij to Eq. ©, 
following the spirit of |32|, • Then, we can notice that 
(oijdxjUj) appears as a production term due to turbu- 
lence for the mean polymer elastic energy. Now, from 
the definition of polymer elastic energy Eq. ((6]), it is 
evident that (E p ) oc ln(/(C fcfc )). So, the FENE-P dumb- 
bells are stretched when —epS^/u^ < in Fig. [P2b and 
then clastic energy is stored on polymers, absorbing tur- 
bulent kinetic energy from the flow. Hence, a mechanism 
of drag reduction can be proposed based on the polymers 
stretching or in other words the behaviour of viscoelastic 
dissipation as a function of the distance y from the wall. 

According to Figs. [T2b and [13] there are three main 
regions in the profiles of viscoelastic dissipation 




0<y/S <di(We C) Re c ) 
<5i(We c , Re c ) < y/S < S 2 (We c , Re c ) 
S 2 (We c ,-Rje c )<y/S< 1. 

(22) 

The first region is at the proximity of the wall, where 
polymers unravel because of the high mean shear, con- 
sistent with other studies [H, [56l - [58l |. storing elas- 
tic potential energy. The range of this region has a 
weak dependence on Weissenberg and Reynolds number 
with its upper bound being within the viscous sublayer 
<5i(We c , Re c ) < 0.05 for all We c and Re c cases considered. 

The second region is the most interesting since poly- 
mers release energy back to the flow, contracting towards 
their equilibrium length, as they are convected away from 
the wall by the near- wall vortical motions. The mani- 
festation of turbulence production by polymers can be 
interpreted in terms of the correlation of the polymers 
with the local fluctuating strain rates and their persis- 
tence in this region. In particular, -(cr^^.u'J reveals 
that —(a'udxU 1 ) as well as —(a' 13 d z u') are responsible for 



the contraction of the dumbbells and consequently for 
the release of the stored elastic energy, since they are 
positively correlated in this region away from the wall 
(see Figs. ITSa and fT5b) . This region exists in an in- 
termediate y/S range, whose upper bound <52(We c ,Re c ) 
is strongly dependent on We c and less on Re c values. 
As drag reduction amplifies for larger polymer relaxation 
time scales this positive region expands to a wider y/5 
range, which dominates the nature of ~ep8 v /u^. at MDR 
(see Fig. [Hfc). 

Finally, polymers transported away from the wall get 
also negatively correlated with the persistent fluctuating 
strain rates (see Fig. ITSb ) and are extended in the region 
(52(We c ,Re c ) < y/S < 1, which is a sink for turbulent 
kinetic energy, prevailing the LDR flows. However, this 
region is diminished for HDR and MDR flows (see Fig. 
112b ) due to the interplay between the productive and 
dissipative inherent features of ep , which mainly depend 
on the polymer relaxation time scale and the existence of 
intense velocity fluctuations that are able to stretch the 
polymer molecules. 

The phenomenology of the proposed mechanism shares 
many similarities with various conceptual models of ear- 
lier works @. In this study, the basic idea is that the 
transport of the elastic potential energy, stored by poly- 
mers near the wall, is mainly associated with the polymer 
relaxation time scale. The latter determines the distribu- 
tion of energy away from the wall and as a consequence 
the near-wall turbulence dynamics weaken. Up to this 
point, the mechanism agrees with the interpretation of 
Min et al. [56| . which is essentially confirmed by the 
present illustrative computations. However, the novelty 
here is that this mechanism is valid for higher We c values 
and levels of percentage DR in contrast to Min et al. [37j , 



who claim that it is not valid for HDR/MDR flows bas- 
ing their arguments on their debatable numerical results 
(see also section |HIC| . 

In addition, the refinement of the proposed con- 
ceptual mechanism resides on the reduction of ep to 
{oijdxjUi) and even more on the correlations (a'udxu') 
and {a' 13 d z u'), which are responsible for the turbulence 
production by polymer coils. The existence of a third 
dissipative region away from the wall is also emphasised 
in this mechanism, where polymers, after their contrac- 
tion, are now stretched by the intense fluctuating veloc- 
ity field. This outer region dominates the viscoelastic 
dissipative dynamics of the LDR regime and diminishes 
asymptotically as We c increases but it never disappears. 
Ultimately, this picture along with the anisotropy intro- 
duced into the components of turbulent kinetic energy, 
i.e. E = \{u' 2 + v' 2 + w' 2 ), comprise the drag reduction 
mechanism deduced in this study. 



VIII. CONCLUSION 

This paper is devoted to the polymer dynamics in vis- 
coelastic turbulent channel flow and their effects on the 
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flow, reproducing turbulent drag reduction by DNS us- 
ing a state-of-the-art numerical scheme in wall-bounded 
flows to solve the FENE-P model. The potential of this 
methodology to capture the strong polymer-turbulence 
dynamical interactions allowed /3 values to remain high, 
more representative of dilute polymer solutions used in 
experiments. Even then, higher percentage DR values are 
obtained for given We c than previous numerical studies. 

The effects of L p and Re c on the results support the 
claims for non-universality of the dynamics for interme- 
diate levels of DR between the von Karman and the 
MDR law. The universal MDR asymptote, on the other 
hand, is reached in this study under the combination of 
high polymer extensibility L p with high enough elasticity 
given by large values of We c at a given moderate Re c . 

The experimentally observed distinct differences in the 
statistical trends of the turbulent velocity field, particu- 
larly for u'_i_ (sec Fig. , are clearly identified with 
the current numerical approach in comparison with other 
simulations, most of which do not even approach such a 
characteristic trend. Overall, the peaks of the statistical 
profiles of velocity and vorticity fluctuations shift away 
from the wall as DR increases, in agreement with other 
experimental and numerical studies, indicating the thick- 
ening of the buffer layer. At the same time, v(3-^{u) in- 
creases towards the centre of the channel for higher We c , 
denoting the importance of viscosity away from the wall 
at these moderate Reynolds number DNS. 

Lumley's phenomenology Lumley on the manifes- 
tation of drag reduction is based on the conjecture of 
coil-stretch transition, i.e. exponential full uncoiling of 
polymer molecules, for the build-up of intrinsic viscos- 
ity. However, our numerical results illustrate that the 
onset of drag reduction and even the MDR asymptotic 
state can be reached while (Ckk) <C L p with L p large 
enough, in agreement with the initial claim by Tabor 
and de Gennes Tabor and de Gennes [l(| that even high 
space-time strain rate fluctuations near the wall can only 
partially stretch polymer coils. We also showed that the 
percentage polymer extension is less but the actual ex- 
tension is more for larger L p , amplifying DR. Thus, large 
polymer coils that do not reach their critical full exten- 
sibility should be of interest to experimental investiga- 
tions on scission degradation of polymer chains and drag 
reduction effectiveness. Such macromolecules would be 
less vulnerable to rupture avoiding the loss of the drag 
reduction effect. Besides, they should be able to stretch 
substantially to make a stronger impact on turbulent ac- 
tivity and consequently enhance percentage drag reduc- 
tion. 

The analysis of the conformation tensor field provides 
great insight into the polymer dynamics and their influ- 
ence on the flow. The dominant anisotropic behaviour 
of the mean conformation tensor, i.e. (Cn) 3> (C12) — 
(C33) > (C22}, due to the mean shear in viscoelastic tur- 
bulent channel flow, influences the anisotropy of the fluc- 
tuating flow field. The anisotropy in the HDR and MDR 
regimes is depicted at the small scales of our DNS outside 



the buffer layer and towards the centre of the channel by 

Different asymptotic rates of convergence are observed 
for the conformation tensor components towards the limit 
of infinite Weisscnbcrg number demonstrating the com- 
plex polymer dynamics even in this simplified dumbbell 
model. In the limit Wes — > 00 polymers are consid- 
ered stiff, i.e. Cij — > (Cij), mostly in the main di- 
rections of elongation and the correlations of the fluc- 
tuating conformation tensor and velocity fields in the 
other directions are assumed to remain minimal at this 
limit. Therefore, (a u ) = A 1 ^2{C 12 ) ^{u) and (a 12 ) = 

^2^- (C22) -^j (u) , with A\ -»■ 1 and A 2 ->• 1 in a region 
somewhere between the wall and the centre of the channel 
in that limit. Our numerical results show that A\ — > 1 
in such a region but not A2. A2 on the other hand is 
about contant in the range 0.2 < y/S < 0.7 and shows a 
tendency towards 1 as Wes increases. 

The following theoretical view was stated in this pa- 
per with regards to the controversy over the existence or 
not of Reynolds shear stress at the MDR limit, which 
is of fundamental importance to the dynamics of tur- 
bulence production at this limit. It is conjectured that 
at the MDR limit {012) is negligible. This was based 
on the idea mentioned above about the stiffness of poly- 
mers at Wes - ► 00 plus the assumption that (C22) be- 
comes negligible at the same limit. Then, it is supposed 
that this behaviour is also valid under both the infinite 
Weisscnbcrg and Reynolds number limits by taking care- 
fully these limits, so that we go along the universal MDR 
asymptotic line. Hence, one might expect an interme- 
diate region 5 V <C y <C S of approximately constant 
Reynolds shear stress, i.e. —{u'v')/v%. — > 1, implied by 
the balance of shear stresses when taking the limits of 
y/S — > and y/S v — > 00 with the reasonable assumption 
that i>j3j^(u) — > for y ^> S„. In summary, the classical 
turbulence generation by —{u'v') seems to survive at the 
MDR limit, based on the above assumptions. 

Polymer-turbulence dynamical interactions were ex- 
pressed through viscoelastic dissipation sp = {o'[jd Xj u! i }, 
which can cither dissipate or produce turbulent kinetic 
energy. For HDR and MDR flows, J ep dy becomes vital 
in the flow dynamics in proportion to J V dy and J En dy 
due to the vast inhibition of Reynolds shear stress and 
fluctuating strain rates, respectively. In particular, a 
different view of the near-wall dissipation dynamics is 
shown for HDR/MDR flows, with the maximum dissipa- 
tion arising away from the wall. It is intriguing to note 
that £p follows a transitional pattern from LDR to HDR 
regime (see Fig. [T2b ) similar to u' + (see Fig. [5^,) and 
(C22) (see Fig. [5^). This characteristic behaviour is also 
reproduced on average in J epdy, where its dissipative 
feature enhances in the LDR regime but attenuates for 
HDR/MDR flows, with the productive nature dominat- 
ing for high percentage drag reduction. Thus, polymers 
get somehow involved in the production dynamics of tur- 
bulent kinetic energy. 
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In view of the current viscoelastic DNS the following 
conceptual picture of drag reduction is deduced, which 
is an extension to and refinement of the mechanism pro- 
posed by Min et al. [ID]. Polymers in the near- wall re- 
gion extract energy from the flow due to the uncoiling 
caused by the mean shear and release some portion of this 
stored elastic energy back to the flow by contracting as 
they move away from the wall. This transport of energy 
depends on Weisscnbcrg number which determines the 
distribution of energy away from the wall. Ultimately, 
this process undermines the dynamics of near-wall tur- 
bulence. Note that polymers also unravel due to velocity 
fluctuations, as they move towards the core region of the 
flow, extracting again energy from the flow. This mech- 
anism appears to be valid for all drag reduction regimes 
with the dissipative and productive elements of viscoelas- 
tic dissipation competing in different parts of the flow for 
different levels of DR. We also observe that correlation 
WijdxjV-'i) is able to resemble the dynamics of ep and 
specifically that {a'i r d x u') and (cr' 13 d z u') are the correla- 
tions responsible for the production of turbulent kinetic 
energy by polymers. 

So far, in the limited context of the FENE-P model 
and at moderate Reynolds number DNS, the proposed 
phenomenology agrees with the majority of experimental 
and numerical data, where dampening of near- wall turbu- 
lence has long been speculated with various analyses and 
interpretations. Here, however, the transfer of energy 
from the flow to the polymers, its redistribution by the 
latter in the flow field and the prevalence of anisotropy 
over the components of E = ^(\u'\ 2 ) in the three Carte- 
sian directions is suggested as a possible cause of drag 
reduction. 
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Appendix: Computational method 

1. Numerical method for the FENE-P model 

The numerical scheme adapted here for non-periodic 
boundary conditions was initially developed by Vaithi- 
anathan et al. 



Vaithianathan et al. [21| for periodic do- 
mains. The main idea behind the high-resolution cen- 
tral schemes employed here is the use of higher-order 
reconstructions, which enable the decrease of numerical 
dissipation so as to achieve higher resolution of shocks. 
In essence, they employ more precise information of 
the local propagation speeds. A key advantage of cen- 
tral schemes is that one avoids the intricate and time- 
consuming characteristic decompositions based on ap- 



proximate Ricmann solvers [20|. This is because these 
particular schemes realise the approximate solution in 
terms of its cell averages integrated over the Riemann 
fan (see Fig. [T5| . 
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FIG. 16: Central differencing approach - staggered integra- 
tion over a local Riemann fan denoted by the dashed-double 
dotted lines. 

Considering the discretisation of the convection term of 
the FENE-P model only in the x-direction, using the re- 
construction illustrated in Fig. 1161 the following second- 
order discretisation is obtained 
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Similarly, Eqs. (|A.2[) - (| A.4|) can be rewritten for 
H™_ 1 , 2 j k - The appropriate choice of the derivative dis- 
cretisation in Eq. (|A.4|) limits the slope so that the SPD 
property for C is satisfied. The SPD criterion for this 
choice is that all the eigenvalues of the conformation ten- 
sor should be positive, viz. A; > and subsequently all 
its invariants should be positive for at least one of the dis- 
cretisations. Note that just det(C) > 0, is not sufficient 
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to guarantee the SPD property for the tensor [59j . In case 
none of the options in Eq. (|A.4|) satisfy the criterion, then 
the derivative is set to zero reducing the scheme to first 
order locally in space. The proof for C being SPD us- 
ing this numerical scheme can be found in Vaithianathan 
ct al. [2l| . The eigenvalues of the conformation tensor in 
this implementation are computed using Cardano's an- 
alytical solution [6(| for the cubic characteristic polyno- 
mial avoiding any complicated and time-consuming linear 
algebra matrix decompositions and inversions for just a 
3x3 matrix. Ultimately, the advantage of this slope- 
limiter based method is that it adjusts in the vicinity 
of discontinuities so that the bounds on the eigenvalues 
cannot be violated, eliminating the instabilities that can 
arise in these types of calculations, without introducing 
a global stress diffusivity. 

The complicated nature of the slope-limiting procedure 
raises difficulties in the case of wall boundaries for a chan- 
nel flow computation, leading to loss of symmetry in the 
results. This had not been encountered by Vaithianathan 
et al. Vaithianathan et al. [61( , since they only considered 
periodic boundary conditions. So, the implementation of 
the numerical method near the walls of the channel was 
modified for this study considering ghost nodes beyond 
the wall boundaries to keep the original formulation un- 
altered, preserving in that way the second-order accuracy 
at the boundaries. The values at the ghost nodes were 
linearly extrapolated from the interior solution (20j . i.e. 
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The time advancement is done simply using the for- 
ward Euler update, treating implicitly the third, the forth 
term on the left hand side and the right hand side of Eq. 
([5]) due to the potential finite extensibility of the poly- 
mer. Hence, the fully discretised form of the FENE-P 
model is 



so that the convection term and the explicit term coming 
from the time derivative can be assembled in a convex 
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where all coefficients s/ > satisfy Yli=i s i = 1> with 
C* being SPD if the matrices C\ are SPD, ensuring the 
finite extensibility of the dumbbell, i.e. the trace of the 
conformation tensor is bounded trC = A1+A2+A3 < L 2 p 
[2l| . The following CFL condition needs to be satisfied 
for the coefficients s/ to be non-negative 
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and it also determines the time step At. Note that this 
CFL condition is more strict than the one for compact 
finite differences (2(| used for Newtonian turbulence com- 
putations. 

The numerical solution of Eq. (|A.6[) is carried out by 
first rewriting it in a Sylvester- Lyapunov form (62j . sep- 
arating the implicit and explicit terms, i.e. 



• I)x = b (A.10) 
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andB = C* + ^ J are 3 x 3 matrices, (I <g) A T + A T <g> I) 
is a 9 x 9 matrix and x = vec(X), b = vcc(B) are 
9x1 vectors. The formula on the right hand side of 
Eq. (|A.10[) can be reduced from 9x9toa6x6 system 
of equations considering the symmetry of the conforma- 
tion tensor. Note that Eq. (|A.10[) is non-linear and can 
now be solved using conventional methods. In this study, 
the Newton-Raphson method for non-linear systems was 
applied using the LU decomposition for the inversion of 
the Jacobian KM K33f . 



(-in-\-l (-in 



with 



At 



\Hi+l/2,j,k Hi-l/2,j,k) 



AaT l+1/2J ' fe «- 1 / 2 »J. 

-^J-( H uj + l/2,k ~ H i.j-l/2,k) 





i-k 1 j 


(A.6) 


Q^i+l/2,j,k + 


h 

-l/2J,fc 




+ C i,j+l/2,k + C i 


h 

,j-l/2,fc 




+ C i,j,k+l/2 + C i 


j,k~l/2> 


(A.7) 



Time advancement 



After obtaining the new update of the conformation 
tensor C"^, the two-step, i.e. three time-level, second- 
order Adams-Bashforth/Trapezoidal scheme is used for 
the time integration of Eqs. through the following 
projection method [64| 

11* — n n 1 1 

= _(3F" - F n ) + -(P; 1+1 + P n ) (A.ll) 



At 



-Vp 



~n+l 



(A.12) 



where 



F = -- [V(u®u) + (u- V)«] + — Am (A.13) 
2 Re c 
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and 

p =r^ v -(zF1^ c -') (A14 » 

with 

r+1 = h[ n+lpdt - (A - 15) 

The incomprcssibility condition V • u n+1 = is verified 
by solving the Poisson equation 

V ■ 11* 

V • Vp n+1 = -^r— (A.16) 



which is done in Fourier space |25[ . It is well known that 
these multistep methods are not self-starting and require 
a single-step method to provide the first time level [HI, 
[65j . In this study, explicit Euler was chosen for just the 
first iteration of these computations, viz. u n = u n ~ 1 + 
AtF n -\ 
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